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Abstract 


2 . Analysis 


An LU implicit multigrid algorithm is developed 
to calculate three-dimensional compressible viscous flows. 
This scheme solves the full three-dimensional Reynolds- 
Averaged Navier-Stokes equation with a two equation k~e 
model of turbulence. The flow equations are integrated 
by an efficient, diagonally inverted, LU implicit multigrid 
scheme while the k - e equations are solved, uncoupled 
from the flow equations, by a block LU implicit algorithm. 
The flow equations are solved within the framework of the 
multigrid method using a four grid level W-cycle, while the 
lc — e equations are iterated only on the finest grid. This 
treatment of the Reynolds- Averaged Navier-Stokes equa- 
tions proves to be an efficient method for calculating three- 
dimensional compressible viscous flows. 


The three-dimensional Reynolds- Averaged Navier- 
Stokes equations, with the Boussinesq eddy-viscosity for- 
mulation, are written in divergence form and then trans- 
formed from the Cartesian coordinate system (x, y, x) to 
the generalized system (&»7,r). The resulting equations 
can be written: 

dW dF dG a# _ 

dt dr] d$ 

dF v dG v dH v (1) 

d£ + dr) d<; 

with 


1. Introduction 

As the availability of larger and more powerful comput- 
ers increases, so will the attention being directed towards 
computing complex three-dimensional compressible viscous 
flows. At present the most sucessful approach has been to 
solve the Thin-Layer Navier-Stokes equations with an alge- 
braic turbulence model. Most notable has been the work of 
Pulliam and Steger 1 and their highly sucessful ARC codes. 

Attempts to produce efficient numerical schemes for 
the calculation of compressible viscous flows have led to 
the Diagonalized ADI scheme of Chaussee and Pulliam 2 
and more recently to the multigrid Runge-Kutta scheme of 
Jayaram and Jameson. 3 As emphasis increasingly focuses 
upon calculating complex three-dimensional flows, the need 
to numerically solve the Reynolds- Averaged Navier-Stokes 
equations efficiently is correspondingly apparent. 

The present work deals with the development of an ef- 
ficient LU implicit multigrid scheme for the numerical solu- 
tion of the Reynolds- Averaged Navier-Stokes equations and 
the two equation k — e turbulence model. The Reynolds- 
Averaged Navier-Stokes equations, which for mathematical 
closure require the modelling of the Reynolds stress tensor, 
are solved by the Diagonally Inverted LU Implicit Multi- 
grid scheme developed by Yokota, Caughey, and Chima. 4 
A Boussinesq eddy-viscosity formulation is used to model 
the Reynolds stress term, where the turbulent viscosities 
are calculated from a standard high Reynolds number k — s 
model. 
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for which the production rate of the turbulence kinetic en- 
ergy is defined as 
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The turbulence modelling constants are chosen to be the 
standard Launder and Spalding 9 values of 


= 0.09 ti = 1.44 t 2 = 1.92 


a k = 1.0 a t = 1.3 

and were never altered during the course of this work. 
(Isotropy is assumed since is a scalar constant rather 
than a non-constant vector quantity.) The transformed 
Navier-Stokes and k — e equations are discretized by a finite 
volume formulation that approximates the spatial differ- 
ences as a net flux across the faces of each mesh cell. Global 


conservation and the admission of possible uniform flow so- 
lutions are insured by evaluating both the inviscid and vis- 
cous flux vectors on the faces of the boundary-conforming 
mesh cells. This procedure requires that the flow and k — s 
variables, viscous stresses, and flux-embedded geometric 
quantities be defined on the faces of the mesh cells during 
the flux evaluations, although it is the cell-averaged flow 
and k — e variables that are calculated during the time and 
spatial marching. The viscous stresses and geometric quan- 
tities are evaluated directly on the cell faces, while the flow 
and k — c variables are averaged over values found in adja- 
cent cells. The unsteady equations can be discretized into 
an implicit approximation that, when written in a linearized 
delta form, produces a numerical scheme whose steady state 
solutions are independent of the time step size used in the 
time marching. The computational effort required to con- 
struct this implicit approximation is kept at a minimum 
by treating only the inviscid fluxes. The delta form, pro- 
duced by linearizing the changes in the inviscid flux vectors 
through a Taylor series expansion about a time level n, can 
be written for the flow equations as follows: 

[/ + n>M{6tA + 6„B + 6 t C)]A&$ k = 

-A t{S ( (F - F v ) + 6„(G- G v ) + $ ( {H - #„))”* (4) 

where 

a w n = w n + l - w n 

and At is the time step size; 0 < < 1 is a parameter 

governing the degree of implicitness; 6 and 5 are cell- and 
face-centered central differences; / is the identity matrix; 
and A, B, and C are the inviscid flux Jacobian matrices 
relative to the vectors F, G, and ff. The implicit form of 
the k — e equations can be written as follows: 

[7 + fi x At(6^Af ct + 8 n B k€ + 6<C k€ - £*«)] = 

— At(6^P kt -f S n G kt + 8<.H kt — S kt )? jk (5) 

and 


A W? € = - W? t 

where A kti B ktt and C kt are the inviscid flux Jacobian 
matrices relative to the inviscid terms found in the vectors 
F kty G kty and H kt \ and E kt is the Jacobian matrix relative 
to the source vector S ks . 

3. Artificial Dissipation 

The finite volume formulation reduces to a central dif- 
ference approximation on a uniform grid and therefore re- 
quires the addition of explicit artificial dissipation terms 
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2) Upper Sweep 


5, Diagonally Inverted LU Factorization 


[/ + *At(# + «+C a )] AW"* - A Y? jk (11) 

A similar explicit sweep, but in directions opposite to those 
taken in the lower sweep, is needed to solve this system of 
equations. The resulting flow field corrections are then used 
to update the flow field. 


Thp LU factorization of the k - e equations produces 
simple 2 x 2 matrix systems which can be inverted alge- 
braically, while the matrix systems associated with the 
flow equations are diagonally inverted using the similarity 
transformation 4 that produces the following diagonal ma- 
trix: 




+ 5 * 


Q~ l {A + B + C)Q = A 


The LU factorization requires the solution of two block 
triangular operators each of which, through back substitu- 
tion, can be reduced to simple 5x5 matrix systems at every 
mesh cell. These reduced systems can be written: 

1) Lower Sweep 
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A similar analysis for the k — e equations can be used 
to write the following 2 x 2 matrix equations: 

1 ) Lower Sweep 
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which has elements 

An — A 22 = A 33 = U 4 V 4 W 


A „ = U + V + W- c\ftf + £+lZ 

A55 = U + V + W + cy/ll + q + ti 
= 0 when i ^ j 

where 

= U 2 4 V 2 4 u; 2 
*1 = £* + 'lx + fx 
h = Cy 4 TJy 4 £y 
h = 6 4 Vz 4 Cz 

and c is the local speed of sound. This similarity transfor- 
mation, for a local time step defined as: 

Cn 

(|A| + |B| + |C|) 

where Cn is the Courant number used in the time marching, 
can be used to transform the lower and upper sweeps into 
scalar equations with the following vector components (m = 
1,...,5) 
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7. Steady State Calculations 

The calculation of steady state solutions is made more 
efficient by using local time-stepping and the multigrid 
method. 

Local time-stepping is used to optimize the time step 
throughout the flow field. A locally varying time step size, 
based on a constant Courant number, is used to create a 
warped time integration that can accelerate the calculation 
to a steady state without affecting the steady state solution. 
These time steps are defined identically for both the flow 
and k — e equations, although the value of the Courant 
numbers used in their evaluation may vary. 

The multigrid method is incorporated into the Diago- 
nally Inverted LU scheme to accelerate the removal of low 
frequency errors from the flow solution and thus increase 
the efficiency of the time-marching procedure. Following 
the work of Jameson, 16 the flow solver is used to smooth 
out high frequency errors resolvable on any current grid 
level (h), while the multigrid method is used to eliminate 
low frequency errors through a sequence of flow calculations 
on coarser grids (2 A, 4/i, 8/i, . . .) . The multigrid sequencing 
used is the four-level W- cycle shown below 



Coarse grid boundary conditions are identical to those used 
on the fine grid with the exception of the inflow/outflow 
conditions which are updated only on the fine grids. The 
Reynolds- Averaged Navier-Stokes equations are solved on 
the finest grid while only the Euler fluxes are evaluated on 
the coarse grids. Coarse grid residuals are kept smooth by 
adding only a constant coefficient second difference artifi- 
cial dissipation term (the nonlinear blended terms are used 
only on the fine grid). This treatment attempts to limit 
the amount of high frequency errors reintroduced into the 
flow field by the upward interpolation of the coarse grid 
corrections. The flow solver, in this case the Diagonally In- 
verted LU scheme, is invoked only once on each grid level 
and only before transfering the flow field to the next coarser 
grid. The multigrid cycle defined above requires approxi- 
mately 1.32 work units of computational effort, where work 


units are normalized by a single Navier-Stokes calculation 
on the finest grid. The k - e equations are solved only on 
the finest grid and are not accelerated with the mtlltignd 
method. Preliminary attempts at multigridding the k - £ 
equations proved to be less efficient than simply solving 
these equations only on the finest grid. 

8. Results 

Numerical results are presented to illustrate the Di- 
agonally Inverted LU scheme’s ability to calculate three- 
dimensional compressible viscous flows and the convergence 
acceleration produced by the multigrid method. 

Turbomachinery calculations were performed on H- 
type grids consisting of 96x24x24 mesh cells in the through- 
flow, blade- to- blade, and radial directions, respectively. 
These grids were generated using a modified version of the 
GRAPE code 17 originally developed by Sorenson. 18 All cal- 
culations were performed on a CRAY X-MP, where a calcu- 
lation consisting of 211 work units required approximately 

1.5 million words of memory and 40 minutes of CPU time. 

The test case used to evaluate the LU scheme is the 
Annular Cascade designed and extensively tested at NASA 
Lewis. 19,20 The computational grid shown in Figure 1 is 
based on the full annular ring of 36 core turbine stator 
vanes. The geometry is a 38.10 mm high untwisted blade 
of constant profile with an axial chord of 38.23 mm. The 
stator has a tip diameter of 508 mm and a 0.85 hub-to-tip 
radius ratio. Mesh cells found immediately adjacent to solid 
walls are centered at distances 0.002 of an axial chord away, 
which correspond to a value of & 60 for the following 
flow calculations. 

Experimental test conditions of ambient axial inflow 
and a 0.65 hub-static to inlet-total pressure ratio produce 
a flow field with mean radius inlet and exit critical velocity 
ratios of 0.231 and 0.778 respectively. To match the up- 
stream flow conditions (an inflow Mach number of 0.211), 
the nonrotating (ROT = 0) calculations were run with a 
0.665 hub-static to inlet-total pressure ratio (PR = 0.665). 

Figure 2 shows the convergence history of the Annu- 
lar Cascade calculation (with multigrid) where a drop of 

3.5 orders of magnitude in the flow field error was pro- 
duced within 200 work units. The residual was reduced 
from an initial RE Si = 0.169Z10 1 to a final value of 
RES2 = 0.572xl0~ 3 (0.9628 average rate of convergence). 
A Courant number of Cn = 6 was used in the solution of the 
flow equations while a Cn = 4 was used in the k — e equa- 
tions (no attempt to optimize these numbers was made). 

The resulting flow field is fully subsonic and is com- 
pared with experimental data at three spanwise positions. 
Figures 3, 4, and 5 compare the calculated blade surface 
static pressure distributions (normalized by inlet total pres- 
sure) at 13.3, 50, and 86.7 percent span with the exper- 
imental data produced by Goldman and Seasholtz. 20 The 
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Fig. 1. Computational grid for the Annular Cascade. 
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Fig. 3. Static pressure blade distribution at 13.3 percent 
span. 
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Fig. 2. Convergence history for the multignd calculation. 
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